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ABSTRACT 

In this paper we make a phase dependent study of the effect of the distortion of 
local magnetic field due to confinement of accreted matter in X-ray pulsars on the 
cyclotron spectra emitted from the hotspot . We have numerically solved the Grad- 
Shafranov equation for axisymmetric static MHD equilibria of matter confined at the 
polar cap of neutron stars. From our solution we model the cyclotron spectra that 
will be emitted from the region, using a simple prescription and integrating over the 
entire mound. Radiative transfer through the accretion column overlying the mound 
may significantly modify the spectra in comparison to those presented here. However 
we ignore this in the present paper in order to expose the effects directly attributable 
to the mound itself. We perform a spin phase dependent analysis of the spectra to 
study the effect of the viewing geometry. 
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1 INTRODUCTION 

Neutron stars in high mass X-ray binaries have high mag- 
netic fields (~ 10 12 G) and accrete matter from their com- 
panion stars either via stellar winds or by disc accre- 
tion. Magnetospheric interaction with the accretion flow 
causes the matter to be channelled to t he magnetic 
poles, forming accretion columns (see e.g. Ghosh et al.l 
il 19771), iGhosh fc Lambl l| 19781 ), iKoldoba et all (|2002l ) and 
iRomanova et all (l2003h ). The infalling plasma, with ini- 
tial relativistic infall velocities, passes through an accretion 
shock at a height of a few kilometres from the neutron star 
surface an d then settles down to a gradually slowing su b- 
sonic flow (|Brown fc Bildstedfl998l : ICumming et al.ll200ll ). 

Such X-ray binary systems show characteristic cy- 
clotron resonance scattering features (CRSF) in their spec- 
tra resulting from resonant scattering of radiation by elec- 
trons in the presence of strong magnetic field (for discus- 
sion on theory and observation o f cyclo t ron scattering fea- 
tures se e e.g. | Harding fc Preecel lll987l). lArava fc Harding 



1999), lAraya-Gochez 



Schonherr et al 



2007) and iMihara et alj 1)20071) ). In the immediate post 
shock region the flow velocities are still relativistic (~ 0.16c 
for 7 = 5/3 gas) and the plasma is optically thin to cyclotron 
scattering. As the accreted plasma descends and cools, it 
forms at the base a static mound confined by the magnetic 
field, and becomes optically thick to cyclotron scattering. 
Any distortion of the magnetic field in the mound due to 



pressure from the confined plasma will be reflected in the 
spectra emitted from the boundary of this region. 

The nature and variation of the cyclotron spectra can 
give important clues regarding the properties of the emission 
region. Many systems show variations of line energies of the 
CRSF with the phase of rotation e.g. Vela X-l, Her X-l, 4U 
0115+63, GX 301-2 etc. This can be due to the variation of 
the local magnetic field structure at one or both poles as a 
line of sight moves across the neutron star. Apart from the 
spin phase dependence, the cyclotron spectra are also seen to 
depend on the l uminosity state of the s ystem. Some systems 
like V0332+53 (|Tsvgankov et al.l l2010) show a negative cor- 
relation between l uminosity and cyclot ron line energy while 
some like Her X-l (IStaubert et alj 20071 ) show a positive cor- 
relation. Such dependence of the line energy with change in 
accretion rate suggests change of local geometry or magnetic 
field structure. Some sources (e.g. 4U 1538-52, A 0535+26, 
V 0332+53 etc) show multiple absorption features with an- 
harmonic separation which can b e due to distortion of lo cal 
field from dipolar magnetic field (|Nishimurall2005l . I201l1 ). 

In this paper we examine the effect on the cyclotron 
spectra arising from the distortion of local magnetic field 
caused by the confined plasma. We consider an accreted 
mound in static equilibrium confined by the magnetic field 
at the magnetic pole of a neutron star. We construct the 
equilibrium solution by solving the Grad-Shafranov equa- 
tion. We do not consider the effects of continued accretion 
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in this paper. We model the X-ray emitting hotspot as a 
mound of accreted matter with finite height and no atmo- 
sphere. The Grad-Shafranov equation for the accreted mat- 
ter on the neutro n star poles has been previously solved by 



ter on the neutro n star poles has been previously solved by 
other au thors e.g. Hameurv et al. (1983'). Br own fc, Bildstenl 
ill998h. ILitwin et al.1 teOOlft. iMelatos fc Phinnevl (12003 )7 



Payne fc Melatosl ||2004) . [Payne fc Melatosl (|2007l ) and 
Vigelius fc Melatosl ( 20081 ) whose main aim was to study the 
extent of deformation and stability of the confined mound 
and also to deduce the effects of magnetic screening on the 
dipole moment of a neutron star. In this paper we extend 
this body of work to predict the cyclotron spectra emanating 
from such mounds. 

We adopt a geo m etry similar to t hat used by 



Hameurv et all (ll983T l. iBrown fc Bildstenl <|l998l ) and 
Litwin et al.1 1|2001^ and a numerical a lgori thm similar to the 
one de veloped by Mouschoviasl (11974 ) and lPavne fc Melatosl 



112004 ) (PM04) for solving the Grad-Shafranov equation. 
However our treatment differs from PM04 in several aspects. 
We work in an axisymmetric cylindrical coordinate system 
instead of the spherical coordinate system of PM04. We use 
a polytropic equation of state for the accreted gas instead 
of the isothermal equation of state of PM04. Finally, we 
consider the mound to be strictly confined to the polar cap 
region, while PM04 allowed a significant amount of mass 
loading outside the polar cap. 

We simulate the cyclotron spectra emitted from the ac- 
creted mound and perform a phase resolved analysis of the 
emission. Our main objective in this paper is to perform a 
phase dependent study of the effects of accretion induced 
distortion of the local magnetic field on the emergent spec- 
tra. In this work we do not perform a detailed radiative 
transfer calculation of CRSF. Instead we use a Gaussian 
profile for the cyclotron feature originating from each point 
of the emission region, with the central line energy given by 
the magnetic field strength at that poin t, according to the 
well k nown relativistic formula given bv lSokolov fc Ternovl 
( 1968). We also incorporate the effects of gravitational bend- 
ing of light and finite energy resolution of detectors. We 
generate the resultant spectra by integrating over the en- 
tire mound, taking into account the variation of the field 
strength over the emitting region. 

We structure the paper as follows. In Sec. ([2)) we first 
review the formulation of the Grad-Shafranov equation for 
the static MHD equilibria. We then outline the numerical al- 
gorithm adopted to solve the Grad-Shafranov equation and 
the test cases for verifying the code. In Sec. [3] we discuss 
the nature of the solutions obtained by solving the Grad 
Shafranov equation and discuss the range of parameter space 
within which the valid solution can be obtained. Our results 
are indicative of the onset of MHD instabilities beyond this 
boundary. In Sec. |4]we describe the algorithm used to simu- 
late the spectra from the mound and discuss the results from 
our simulation of the cyclotron absorption features. In Sec. 
[5] we summarise the results obtained from the simulations 
of the spectra and discuss the implications on observations 
of actual sources. The technical details of the geometrical 
construct used to compute the spectra are presented in Ap- 
pendix [B] 



2 STATIC MHD EQUILIBRIA OF ACCRETED 
MATTER ON NEUTRON STAR POLES 

In this work we consider a neutron star binary system where 
the ma gnetosphere cuts off the accretion disc at Alfven 
radius (|Ghosh et al.lll977l ; [Ghosh fc Lamblll978f ) and mat- 
ter is accreted on to the polar cap. We will consider a 
typical slowly spinning neutron star of mass 1.4Mq, ra- 
dius R = 10 km and magneti c field B = 10 12 G (e.g. 
iBhattacharva fc van den Heuvell (|l99ll )V A polar cap of ra- 
dius R p = 1 km will be considered corresponding to the 
footprint of dipole field lines that extend beyond a typical 
Alfven radius of ~ 1000 km. Accreted matter is assumed to 
be confined within the polar cap region. We will consider 
the accreted matter to form a degenerate mound of finite 
height (~ 100 m or less) with a polytropic equation of state. 
We assume that the mound is in a steady state equilibrium 
supported by the magnetic field. We work in a cylindrical ge- 
ometry (r, 9, z) with origin at the base of the polar cap (see 
Appendix iBl) and consider Newtonian gravity with co nstant 
acceleration |Hameurv et al.lll983l ; ILitwin et al1l200ll ). 
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The initial magnetic field (when no accreted matter is 
present) is dipolar. We approximate the dipolar field in the 
region by an uniform field along z : B<j = Bqz. We as- 
sume axisymmetry of the polar cap mound and use the ideal 
MHD equations, which may be cast in the form of the Grad- 
Shafranov equation. We solve this numerically to find the 
field and matter density configuration for the static equilib- 
rium solution of the system. 



2.1 Formulation of Grad-Shafranov equation 

For an axisymmetric system, one may decompose the mag- 
netic field into a poloidal and a toroidal part : 



B = B,, 



B« 



Vip x 



(2) 



For our work we will assume Bg = 0. The function ip(r, z) 
is the flux function which at a fixed r and z is proportional 
to the poloidal flux passing thr ough a circle of radius r (see 
iKulsrudl l|2005l ); lBislanrDl (|l993r ) for more discussion on this) . 
The poloidal components of the magnetic field are 



r oz 



r dr 



(3) 



Using eq. ([2]) we can write the static Euler equation as 



Vp - ps + 



4:Tvr'- 



■Vip = 



(4) 



where A 2 is the Grad-Shafranov operator : A 2 = ) 
■^j. Assuming an adiabatic gas p = k a dp 



8 (1 d 



we sep- 



arate the r and z components 1 Hameurv et al.l 1 19831 ; 
Litwin et al.l l200l | ) bv the m ethod of characteristics (sim- 
ilar to lPavne fc Melatosll2004 



Pz- pg + 



aV 



i> z — ; p r + 



A 2 V 



Ipr = 



(5) 



4 7rr 2 Aur 2 
where the subscripts indicate partial derivatives. Eliminat- 
ing ^3 from ([5]) we get the equation of the integral curve 
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dz 



dr 



Ipz/lpr 



dp 



(6) 



where c 2 = "fp/p is the adiabatic speed of sound. Solving the 
above two equations, we get ip = constant (which means the 
solutions are on constant ip surfaces) and gz + 



(7-l)p 



f(ip)- Here f(ip) is a ip dependent constant of integra- 
tion. Rearranging th e terms we can write the density as 
ijHameurv et alJll983h 



p = A[ZoW 



(7) 



where A — [g(y— l)/(7k a( j)] is a constant. The function 
Zo(ip) is the mound height profile which defines the vertical 
height of the mound for an adiabatic gas expressed in flux 
coordinate (ip) instead of r. The values of p,p and their 
derivatives go to zero smoothly at z — Zo(ip). 

Putting eq. Q in eq. ([4]) we obtain the Grad-Shafranov 
equation (hereafter G-S) for an adiabatic gas. 



A 2 ip 
4nr 2 



-pa 



d_Zo 

dip 



(8) 



The Grad-Shafranov equation is a coupled non-linear ellip- 
tic partial differential equation. We have solved the Grad- 
Shafranov equation numerically following the algorithm out- 
lined in Appendix [X] 



3 RESULTS AND ANALYSIS OF SOLUTIONS 
OF GRAD-SHAFRANOV EQUATION 

We have made several runs with different mound height pro- 
files. The solutions show expected behaviour of matter push- 
ing field lines outwards until the tension in the field lines 
support the gas pressure. The deformation of the field lines 
increase with larger base pressure and density of the mound 
till the solution breaks down after a threshold density (see 
Sec. 13. 3[) . In this section we discuss the solutions from some 
sample runs and the range of parameters for which equilib- 
rium solutions can be found. 



3.1 Modelling the magnetically supported 
accretion mound 

We will assume a hydro gen poor plasma (/j e = 2) being 
confined in the mound |Brown fc Bildstenl Il99ct ). We re- 
strict the analysis to the gaseous state before ions form a 
liquid phase. The electrostatic coupling parameter (F) gives 
a rough estimate whether matter is solid (r > > 1), liquid 
(r ~ 1) or gas (r < 1) (e.g. lLitwin et all l|200ll )) 



r 



ZV ,47rn. 



1/3 



k B T y 3 
Z 2 



1.1 



A 1 / 3 



10 8 g cm - 



1/3 



10 8 iY 
T 



(9) 



Hence mounds of base densities < 10 s g cm -3 are con- 
sidered for this work. To determine the appropriate form 
of equation of state, we first check if the plasma is non- 
relativistic (7 = |) or relativistic (7 = |) by evaluating the 
adiabatic index (7 = from the expression of fer mionic 

pressure for degenerate electron gas (|Chandrasekharlll967t l. 
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Figure 1. Plot of adiabatic index (7) vs density (g cm" 3 ). For 
lower densities 7 asymptotically converges to 1.6667, which is 
the value obtained in the non-relativistic approximation. For p > 
10 8 g cm -3 7 converges to 1.333 which is the value obtained in 
the ultra-relativistic approximation. 

As shown in Fig. [1] at p ~ 10 7 g cm -3 we have 7 > 1.4. For 
densities lower than this 7 rises and reaches 1.667 asymp- 
totically. Since, for p < 10 7 g cm -3 , 7 is significantly higher 
than 1.33 we model the accreted matter in the mound as a 
degenerate non-relativistic electron gas with a thermal pro- 
ton background whose pressure is negligible compared to 
electron degeneracy pressure. The equation of state for such 
a system is (values quoted are in cgs) 

5/3 



P = [(30 



2 ,.2/3 



K 2 
5m e 



p e m p 



3.122 x 10 1 p° 



(10) 



The plasma is dominated by degeneracy pressure if 
< 1, where Tp is the Fermi temperature : 



m e c 



1-1] 



X F = ^- = — («^ 



-) 1 ^ 3 p 1//3 , pf being the Fermi mo- 
mentum. Obj3erved_2^a2MMnti^^ X-ray 
binaries (|Coburn et al.ll2002l ; iBecker fc Wolfi1l2007h indicate 
that photospheric temperatures in the hotspot are in the 
range T ~ 5 — 10 keV. The Fermi temperature falls below 
10 keV only for a thin layer (~ 0.0LZ o (V')> from «!■ 
at the top and a major fraction of the mound has Fermi 
temperature larger than lOOkeV. Comparison to the tem- 
perature profiles o b tained including heat transfer effects by 
iBrown fc Bildstenl (|l998T ). shows that the temperature in- 
side the mound remains lower than Fermi temperature at 
greater depths. Thus modelling the confined accreted mat- 
ter as a degenerate mound is appropriate. 

We have solved the G-S equati on for differen t form s 
of the mound heigh t profile (e.g. iLitwin et al.l l|200lf ); 
iHameurv et"aH (|l983h ) 



Z C {1 



(f) 2 ) 

Ipp 

Z (iP) = Z c e X p(-2p) 

ip P 



Z c (l 



(11) 

(12) 
(13) 



The mound height profile depends on the mass loading func- 
tion at the accretion disc and redistribution of matter in the 
shock region for which at present there is no clear knowl- 
edge. We resort to evaluating the density by specifying the 
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Density mound on polar cap 
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0.4 

V 

x 

0.2 

0.0 
0. 



0.2 0.4 0.6 0.8 
Radius in Km 



Figure 2. The shape of the accretion mound plotted along r and 
z axis in equal scale to show the real aspect ratio. The mound is 
like a thin flat layer on the pole. 
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Field line plot for Z c = 55m 
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Field line plot for Z c = 70m 
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Figure 3. Solution for Z c = 55 m (top fig) and Z c = 70m 
(bottom fig), cases (a) and (b) in the text. Solid lines are field 
lines from G-S solution. The dash-dotted line represents the top 
of the mound. 



mound height as a simple function of ip, subject to the con- 
straint : p — > as r — > R p , so that the mound is confined 
within the polar cap. In most of the analysis we have used 
eq. (|1H) which has relatively shallow gradients and helps in 
speeding up the numerical convergence. 



3.2 Solutions from the GS-solver 

We discuss here results from two sample runs 
Case (a) : Z c = 55 m 
Case (b) : Z c = 70 m 

with the mound height profile specified by eq. (|f f } and 



Zc=75m 



s t 

g 0.3 -\ 




Zc=76m 



Iteration steps 



Figure 4. Mean ip as a function of iteration steps for different 
mound heights above the stability threshold. The mean ip is seen 
to oscillate between multiple states. Beyond a certain Z c it passes 
through different states randomly. 



f 12 G. The value Z c was chosen to keep 
g cm- 



magnetic field B ■ 

maximum base densities less than fO 8 g cm -3 (as discussed 



in Sec. EE}- 

For case (a) the total mass of the mound is ~ 9 x 
1O _13 M and maximum base density is ~ 4.7 x 10 g cm . 
For case (b) the total mass of the mound is ~ 2.13 x 
10" 12 M Q and maximum base density is ~ 6.8 x 10 gem . 
The mound is in the shape of a flat thin layer on the surface 
of the star, confined within the polar cap (Fig. E2}. Contours 
of ip from the solution, which represent the magnetic field 
lines (as B ■ V?/i = 0) are plotted in Fig. E3 From the figure 
we see that the field lines are bent to support the pressure 
of the confined matter. The distortion is more in case (b). 
Field lines are pushed outwards from the initial configura- 
tion resulting in bunching of field lines and increase in field 
strength. 



3.3 Valid parameter space for existence of 
solution 



lHameurv et"ai1 (|l983h mention in their work that for the 
configuration of mound height profiles they had considered, 
no solution was found for field lower than a critical value. 
We observe the same for different values of magnetic field 
and different mound height profiles. For a fixed magnetic 
field we find that for the parabolic profile (eq. Ill|l , a stable 
solution exists for only up to a maximum threshold mound 
height (Zmax)- For mounds higher than this, the ip func- 
tion keeps oscillating between multiple states with closed 
magnetic loops during the iteration process and conver- 
gence to an unique solution is not reached. For magnetic 
field ~ 10 12 G the maximum height of a mound for a sta- 
ble solution was found to be Z max ~ 70m. For a mound 
higher than this threshold, Z c — 75m, the ip at intermedi- 
ate steps pass through states similar in nature to states 1 
and 2 as shown in Fig. E5] The mean ip is seen to oscillate 
between two states as in Fig. [4] For higher mounds the 
branches of mean ip bifurcate to multiple states similar to a 
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Field lines for state 1 



Field lines for state 2 




0.8 1.O.0 0.2 
Radius in Km 

Figure 5. The tp function at two intermediate iteration steps 
(79 th and 80 th ) of the GS-solver for a mound of height 75m. 
Closed magnetic loops are seen to form which indicate loss of 
equilibria. At different iteration steps, the ip passes through states 
similar to states 1 and 2 depicted here without reaching conver- 
gence. 



Maximum mound height vs field 

100.0 1 




10 11 
log 10 (B) 

Figure 6. Maximum height of the mound (Z max ) that can be 
supported by different surface field strengths. Mound height pro- 
file of eq. {TTJ is assumed. The dashed line shows a power-law fit 
to the points. Valid solutions can be obtained only for parameters 
below the line, above which G-S code does not converge. 



pitch-fork diagram. Format i on of cl osed magnetic loops, also 
rep orted in|Hameurv et"afl (119831 1 , IPavne fc Melatosf (|2004T i 
and lPavne fc Melatoj (|2007h , appear to indicate loss of equi- 
libria. We have tried different simulations to check for the 
existence of solutions beyond the threshold, e.g. higher reso- 
lution runs or improving the initial guess solution by starting 
from a previously converged solution or increasing the radial 
extent of the grid. However stable solutions were not found 
for heights greater than the threshold Z max . 

Fig. [6] shows the maximum values of Z max (eq. 
for a given field up to which solutions exist. The maximum 
allowed height (i? max ) has a power-law dependence on B 



log 3.676 + 0.4611og 10 (£) 



(14) 



where Z max is in metres and B in Gauss. A similar magnetic 
field to Z max scaling was observed for solution with different 
mound height profiles e.g for eq. (fT2"1) we get log 10 (Z max ) = 
-3.79 + 0.461og 10 (B) and for eq. ([13} we get log 10 (Z max ) = 
—3.61 + 0.451og 10 (B) which indicates that this is a generic 
feature of the Grad-Shafranov equation in the current setup. 

Such a scaling relation can be understood approxi- 
mately by comparing the variation of pressure and magnetic 
field over different length scales which balance each other. 
Lateral variation in pressure over scale R p is balanced by 
tension from curvature in magnetic field which occurs over 



a length scale Z n 
(TTT1) we get 



Hence from eq. ©, eq. (|10|) and eq. 
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Vp ~ B ■ VB 
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(15) 



iLitwin et al.l (1200 ll ) have shown that ballooning instability 
will disrupt the equilibria if A/3 > 7.8_R p /[(7 ~ l)^o (?/>)], 
where /3 is the ratio of plasma pressure to magnetic pressure 



[p/(B 2 /8tt)]. Using p ~ k ad p 5 / 3 (eq. [TOJ and p ~ AZUL (eq. 
[7J| we can write the stability criterion obtained by Litwin as 



l°gio( Z " 



> 



-5.1 + -log 10 



B 



(16) 



which is very close to the observed dependence of Z max and 
B as obtained from our numerical solutions. Thus limit rep- 
resented by eq. (|14p may result from ballooning type pres- 
sure driven instabilities where curvature of magnetic field 
can no longer support the plasma pressure (eq. 1151 eq. I16|l 
and the equilibrium solution cannot be obtained. Hence for 
our analysis of the cyclotron line features in the following 
section, we restrict ourselves to mounds of height less than 
70m for a dipole field of 10 12 G. A detailed study of the 
stability analysis of our solutions will be reported in a fu- 
ture publication (Mukherjee, Bhattacharya and Mignone in 
preparation). 



4 CYCLOTRON RESONANCE SCATTERING 
FEATURES (CRSF) 

The major part of the mound forms an optically thick 
medium with cyclotron line formation taking place in a thin 
layer located at the top surface. The depth of the line form- 
ing region may be estimated as I ~ Zo — z where Zo (the 
mound height profile) is the top height of the mound at a 
given r. From the definition of optical depth and using eq. 
O for the density distribution we find the relation between 
optical depth and thickness of the line forming region as 

T = (17) 



fi e m p 

where m p is proton mass and a is the scattering cross sec- 
tion. For Thomson scattering (a — (7Th), t — 1 occurs at 
a depth of ~ 1.1 mm and for cyclotron resonance scatter- 
ing (a ~ 10 5 CTTh) t ~ 1 occurs at 11.3 fim. Thus cyclotron 
line formation takes place is a thin layer at the top of the 
mound. Variation in the local magnetic field at the top of 
the mound is expected to cause variation in the cyclotron 
spectra. Modelling of the cyclotron resonance scattering fea- 
tures (hereafter CRSF) taking into account the contribution 
of different parts of the mound to the line of sight is pre- 
sented in the following section. 

4.1 Modelling cyclotron spectra 

The emission profile from a point on the mound depends 
on the strength and direction of the magnetic field and the 
angle between the emergent radiation and the local normal 
to the mound surface, which vary with position due to cur- 
vature of the field lines and the shape of mound surface 



© 2011 RAS, MNRAS 000, ??-?? 



6 Mukherjee and Bhattacharya 



respectively. We divide the mound into small area elements 
(AAr^Bj in a (r,0) grid in cylindrical coordinates on the 
polar cap) and find the local magnetic field vector and lo- 
cal normal to the mound for each element, assuming them 
to constant over the area element. The resultant cyclotron 
spectrum is constructed by integrating the weighted contri- 
bution of emission from all parts of the mound towards a 
line of sight (hereafter los) 

F = ^2l(e ctl )AA ri ,e j cosO a i (18) 
id 

I is the normalised intensity at the mound surface and the 
angle 9 a i (eq. IB6[) is the angle between the direction of emis- 
sion at the mound surface and the local normal to the sur- 
face. For integrating the intensity over the mound we have 
chosen a radial grid with resolution equal to or higher than 
the resolution in the radial direction of the G-S solution. 
The cyclotron line energy is given by 

E n = riE c oVT^ (l - | ( 51 ^°y ) sin 2 9 ab ^ (19) 

where n — 1,2,3... is the order of the harmonic, E c o — 
II.6B12 in keV, 8 a t is the angle between the direction of 
emission and local magnetic field (eq. IB5[) and u — r s /r, r 3 
being the Schwarzschild radius. The factor y/l — u gives the 
gravitational redshift of the line. 

Eq. (|19|l is correct to second order in the small param- 
eter _E c o/(511keV), and is adequate for the field strengths 
we consider in this work. In our studies we consider only 
the effect of the fundamental n = 1 line. The accurate de- 
pendence of the width and depth of the CRSF on the angle 
Oab can be obtained by solving radiative transfer equations 
in the mound. This is beyond the scope of the present work 
and will be addressed in a future publication (Kumar et al in 
preparation). For our present work we model the cyclotron 
feature by a Gaussian profile with line centres from eq. (|19p 
and model the dependence on the angle 9 a t, of the line width 
AE/E C and the re lative depth by interpol ating from the re- 
sults presented in ISchonherr et al.l (|2007l ) for the slab 1-0 
geometry. 

We incorporate effects of gravitational light bending 
(eq. IBID foll owin g the approximate formula given by 
iBeloborodovl j2002). For the intrinsic intensity profile we 
use a form I(8 a i) = lo + Ii cos 8 a i, but set Ji = for most 
of our analysis. To simulate the finite energy resolution of 
the detectors, we convolve the spectra with a normalised 
Gaussian with the standard deviation a fraction / of the 
local energy 

w { E,E>) = ^^(dE_^y, a „ fE > ( 20 ) 

Detectors currently used for observations in X-ray astron- 
omy usually have a energy resolution of 10%-20% (/ ~ 
0.1 — 0.2). We carry out the above analysis at different phases 
of rotation of the neutron star to perform a phase dependent 
study of the spectra. The nature of the spectra also depends 
on the relative orientation of the mound (inclination angle 
r) p ) and the los (angle i) with respect to the spin axis of 
the neutron star, which are treated as free parameters (see 
Appendix [B]for details on the geometry). 



Light curve for single mound (Z c = 45m) 
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Figure 7. The light curve of one pole mound for Z c = 45 m, 
rjp = 10° and i = 30° (Appendix [Bj. Insets show the spectra at 
the two rotation phases marked. No significant variation of line 
energy is seen. A detector resolution of 10% was assumed (eq. 

Op. 



4.2 Results : cyclotron spectra from a single 
mound 

We have carried out simulation runs for the cyclotron spec- 
tra using the solutions of the magnetic field obtained from 
the GS-solver. For most of our analysis in this section we use 
the solution with profile eq. (|lip . Although the shape and 
nature of the CRSF will change for different profile functions 
we can draw some general conclusions about the dependence 
of the cyclotron spectra on the local magnetic field. 

We first consider the case of emission from a single 
hotspot at one of the poles. Fig. [7] shows the light curve 
from a mound of height Z c = 45 m at a pole with incli- 
nation rjp — 10° and los i = 30° (See Appendix [B] for 
definitions of i and rjp). The inset plots show the cyclotron 
spectra convolved with a Gaussian with / = 0.1 (eq. I20[) , 
at two rotation phases 32° and 200°. Although the magni- 
tude of the field at the top of the mound varies by ~ 27% 
in this case (Fig. [8]), the line centre of the CRSF shows 
less than 0.2 % change about a mean ~ 9.6 keV. As the 
continuum emission is assumed to be isotropic and uniform, 
and also since gravitational bending redirects the light rays 
in directions well away from straight trajectories, all parts 
of the mound contribute towards a given line of sight at all 
phases. This gives a resultant averaged spectrum with very 
little phase dependence of the CRSF from a single pole. 

For mounds with large distortion of surface magnetic 
field (e.g. for Z c = 70 m, the field at the edge rises to ~ 4.6 
times the original dipole value (see Fig. [8]). We find two 
distinct CRSF fundamentals at two distinct energies (Fig. 
[5]). The feature at a lower energy, for a field ~ 10 12 G, orig- 
inates near the centre of the mound while that at a higher 
energy arises in the high field regions near the periphery 
of the mound (Fig. [H|. We emphasize that this multiple- 
featured spectrum is a result of the variation of the local 
field strength and does not represent multiple harmonics, as 
only n = 1 features have been included in our computation. 
The energy ratio of these features may therefore be arbitrary 
and not follow a harmonic relation. 

The line centres of the individual peaks show little vari- 
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Figure 8. The variation of the strength of the magnetic field 
with radial distance at the top of the mound, for different Z c . The 
plots are offset by an amount (B ff sct ) for clarity. The maximum 
magnetic field at the top is several times the surface dipole field 
due to distortion from pressure of accreted matter. 



Spectra from different mounds 

1.101 i i i | i i i | i i 



1.08 



1.06 



1.02 



1.00 



0.98 





Z c =60m 



20 40 
Energy in KeV 



60 



Figure 9. CRSF from mounds of different heights. The spectrum 
for each mound shows two CRSF fundamentals at two distinct 
energies, corresponding to the undistortcd dipolar field and the 
region with large distortion of the magnetic field respectively (see 
Fig. ®. 
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Figure 10. The light curve of one pole mound with Z c ~ 60 m, 
rjp = 10° and i = 30° (Appendix [Bj. Insets show the spec- 
tra at two the rotation phases marked. Multiple CRSF are seen 
with variation of the relative depth with phase, but no phase de- 
pendence of the line energies. A detector resolution of 10% was 
assumed (eq. I20[l . 



ation with phase but the relative depth of the peaks depend 
on the viewing geometry and the phase angle. Fig. 1 101 shows 
the light curve and spectra (with / = 0.1 in eq. I20[) from 
a mound Z c ~ 60 m with r/ p = 10° and i — 35°, where 
the relative depth of the CRSF vary with phase. The degree 
of their variation depends on the relative orientation of the 
pole and the los, occurring for a certain range of viewing 
angles (25° < i < 40° for a pole at r/ p = 10° in the present 
case). Two distinct CRSF are observed for all mounds of ap- 
preciable magnetic field distortion at the top and all viewing 
geometries, which shows that it is a generic feature of the 
cyclotron lines originating from the mound. 

The cyclotron features are dominated by the field struc- 
tures near the periphery of the mound as these regions have 
a larger emitting area. For G-S solutions with mound height 
profile as in eq. the regions of higher magnetic field 

distortion occur near the edge of the polar cap. Hence we 
observe two CRSF of comparable depths for all mounds with 
profile eq. (|11[) . However, if we use an exponential profile 
(eq. I12|) which falls off sharply with r, regions of larger field 
distortion occur much closer to the axis (Fig. and the 

CRSF is dominated by the undistorted field in the outer re- 
gion. A shallow feature at ~ 18 keV is contributed by the 
higher field regions nearer to the axis. This shows that the 
CRSF is heavily influenced by the geometry and distribution 
of the local field as well as the structure of the mound. 



4.3 Effects of finite energy resolution of detectors 

The finite energy resolution of the detector can make the 
two absorption features indistinguishable if closely placed. In 
Fig. [I2]we show the effect of Gaussian convolution with dif- 
ferent values of / (eq. I20p . For / = 0.2, the two peaks are in- 
distinguishable. For mounds of lower height (e.g. Z c = 45 m) 
the two CRSF become indistinguishable at a much smaller 
/. Thus the finite energy resolution of the detector can often 
mask the effects of the internal magnetic field structure on 
the CRSF. However for mounds of higher height, the two 
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Exponential profile with Z c = 55m 
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Figure 11. Magnetic field line structure (upper panel) for a 
mound with Z c ~ 55m with an exponential mound height profile 
eq. (1121 and the spectrum emitted (lower panel). The spectrum 
is dominated by area elements near the polar cap edge (corre- 
sponding to the region of undistorted field in this case) as they 
span larger areas. The shallow feature at ~ 18 keV corresponds 
to the region of large field distortion closer to the axis. 
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Figure 12. The effect of finite detector resolution modelled by 
convolving the spectra with a Gaussian function. Results for dif- 
ferent values of / (eq. I20|l are shown. The double CRSF nature of 
the spectra disappear for / > 0.2 corresponding to a 20% energy 
resolution of detector. 
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Figure 13. The light curve for two anti-podal mounds, Z c ~ 
45 m and 52 m, with r\ v = 30° and 150° respectively, and i = 30° 
(Appendix |Bj. Insets show the spectra at the two rotation phases 
marked. The CRSF shows 20% change in line energy with phase. 
A detector resolution of 15% was assumed (eq. 1201 . 



CRSF have centres far enough to be distinguishable even 
with / = 0.2. 



4.4 Emission from two anti-podal poles 

Many systems show variation of the cyclotron energy 
with the phase of rotation of the neutron star e.g . 
AE cyc /E cyc ~ 10% for Vela X-l (IKreykenbohm et alfeooa) . 
~ 20% for 4U0115+63 ([Heindl et al.ll2004l; iBaushevfhoog ), 
~ 25% for Her X-l JKlochkov et ail l2008bft . GX 301- 
2 (IKreykenbohm et al.1 |2004|). and~^~ 30% for Cen X-3 
l|Suchv et al.ll2008l ; iBurderi et all 120001 ) . From our simula- 
tions of CRSF for emission from a single hotspot we find 
very little variation of the line energy with rotation phase. 
However if emission is considered from two antipodal host- 
pots at opposite poles with slight difference in mound height 
due to asymmetric accretion, then the line centre of the re- 
sultant CRSF show a stronger phase dependence. Fig. 1131 
shows the light curve for a neutron star with two anti-podal 
poles at rj p ~ 30° and 150°, having mounds of height Z c ~ 
45 m and 52 m respectively, and an observer at inclination 
i = 80°. The simulation is carried out by assuming uni- 
form and isotropic continuum intensity normalised to 1 at 
both poles, which may not be valid in reality due to differ- 
ence in accretion rates at the two poles. However for small 
differences in the mound heights considered here, we ignore 
the differential luminosity effect and draw some general con- 
clusions about the behaviour of the cyclotron line energies. 
The spectra are convolved with a Gaussian function with 
/ = 0.15 (eq. HI). 

The light curve is not sinusoidal unlike the case of a sin- 
gle pole (Fig. and Fig. [10). The line energy of the CRSF 
varies by AE C ~ 20% over a full rotation cycle. The spec- 
trum from the pole nearer to the los dominates the CRSF. 
The observed variation in the line energy of course depends 
on the viewing geometry defined by the angles i and rj p . A 
proper evaluation of the spectra for a two pole case would 
require the knowledge of the accretion rate and local tem- 
perature of the hotspots, but it may be concluded that in 
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Figure 14. The CRSF from a mound of height Z c = 60 m with 
mound height profile eq. 1111 rj p = 10° and i = 70° (see Ap- 
pendix IB[>. The undistorted field produces the small feature at 
~ 8 keV. The distorted field near the polar cap edge is the domi- 
nant contributor to the CRSF, producing a dip at ~ 18 keV. This 
shows that field inferred from CRSF line energy may not depict 
the surface dipole field, but the field distorted due to accreted 
matter. 



the presence of multiple hotspots, the CRSF line energy will 
show a significant phase dependence. 



5 DISCUSSION AND CONCLUSIONS 

We have modelled the structure of the accretion mound 
by solving for the static magnetic equilibria described by 
the Grad-Shafranov equation. For mounds of total accreted 
mass ~ lO _12 A/0 (Sec. I3.2[) there is appreciable distortion 
of the magnetic field. We have found that for a given surface 
field strength, stable solutions to the G-S equation are not 
found for mounds of height higher than a threshold (Sec. 
13. 3|) . Beyond the threshold, field lines with closed loop con- 
figurations are formed, and may indicate the onset of pres- 
sure driven MHD instabilities. So we restrict our analysis to 
mound heights for which equilibrium solution is obtained. 

Using the mound structure obtained from the Grad- 
Shafranov equation, we have simulated the cyclotron spec- 
tra which will originate from its mound surface. We have 
integrated the emission from different parts of the mound to 
find the resultant spectra and have performed a phase de- 
pendent study of the spectra. We have assumed a Gaussian 
model for the CRSF fundamental and have incorporated ef- 
fects of gravitational bending of light and finite energy res- 
olution of the X-ray detectors (Sec. 14. ip . From our analysis 
of the phase dependent spectra (Sec. 14.21 Sec. 14.31 and Sec. 
14.41) we can draw some general conclusions : 

(i) CRSF need not represent the dipole field : The 

CRSF are heavily influenced by the distortion of the local 
field caused by the pressure of the confined accreted mat- 
ter. The line energy of the CRSF may not always represent 
the intrinsic dipole magnetic field. As shown in Fig. 1141 
the small dip at ~ 8 keV is the redshifted cyclotron feature 
due to the surface field of 10 12 G, which may not be ob- 
served in practice in the presence of noise and poor detector 
resolution. The CRSF at ~ 18 keV resulting from the dis- 
torted field will be the dominant feature in the spectrum. 
For mounds of height Z c ~ 70 m the maximum distortion 
in the field can be as large as 4.6 times the dipole value. 



(ii) Multiple, anharmonic CRSF : From our simula- 
tions we observe Multiple CSRF fundamentals from a single 
mound with considerable field distortion (see Fig. O Fig. 
I10[) . A real spectrum will also contain multiple CRSF har- 
monics which will add to the complexity of the spectrum. 
Several HMXBs show multiple CRSF where the line centres 
have anhar monic separation e.g. E c ~ 22 keV, 47 keV for 

23 keV, 51 keV for 
45 keV, 100 keV for 



4U 1538- 52 dRodes-Roca et al 



Vela X-l (iKrevkenbohm et al.l 



2009), 



200: 



A 0535+26 (|Kendziorra et al .ill 994 ; ICaballero et~al.l 1200 
26 keV, 49 keV, 74keV for V 0332+53 (jMakishima et'al 

36 keV, 63 keV for EXO 



1 19901 ; iPottschmidt et al1l2005h, - ■ ^ , . 
2030+375 (|Reig fc Coelll999l ; iKlochkov et all l2008aT ) . More 



on the p roperties of these so urces can be found in the re- 
views by iMihara et al.l (|2007T l and iLutovinov fc Tsygankovl 
(|2008h . 

Anharmonic line spa cing due to i ntrin sic non-dipolar field 
has been discussed by iNishimural l|2005h . Our results show 
that strong non-dipolar structure can be generated in the 
accretion process itself. 

(iii) Effect of detector resolution : The detectability 
of multiple cyclotron features depend on the energy resolu- 
tion of the detector. For detectors with poor energy resolu- 
tion, the multiple fundamental features will be masked as 
shown in Fig. (|12[) . However for mounds with large field 
distortion (Z c ~ 60— 70m) multiple absorption features will 
still be observed with detectors that have energy resolution 
of < 20%. 

(iv) CRSF phase dependence - one pole : We find 
no appreciable spin phase dependence of the line energy of 
a given CRSF from a single mound (Fig. [7} despite there 
being considerable local variation of magnetic field. This 
may be attributed to the fact that all parts of the mound 
contribute towards a given line of sight at any phase due 
to effects of gravitational light bending and assumption of 
isotropic local emission. However we find that the relative 
depth of the multiple CRSF depends on phase and viewing 
geometry as shown in Fig. 1101 

(v) CRSF phase dependence - two pole : Line energy 
variation is however observed if emission from both poles 
with slightly different mound heights are considered (Sec. 
14. 4|) . For emission from two anti-podal mounds of height 
Z c ~ 45 m and ~ 52 m the CRSF line energy varies by 
20% during one spin cycle, similar to what is observed for 
many sources. Thus we conclude that emission from multiple 
accretion mounds will result in strong variation of the line 
energy of the CRSF with spin phase. 

(vi) Dependence on mound structure : CRSF are 
dominated by field structure at the mound periphery due 
to the larger physical area of this region. Different structure 
and density distribution of the mound would cause different 
distributions of the local field. Fig. [9] and Fig. Illl show the 
difference in the CRSF from a mound of the same maximum 
height (Z c = 55 m) but different mound height profiles. The 
strong dependence of the spectra on the structure and size 
of the mound suggests that with variation in the accretion 
rate, the observed line profiles and energies may change, 
contributing to a luminosity d ependence of the spec trum. 

Some sources e.g. Her X-l (Istaubert et al.ll2007h and A 
0535+26 (|Klochkov et al.ll2011I L~ show a positive correlation 
between luminosity and line energy. In our picture, this can 
be attributed to the presence of a strong non-dipolar com- 
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ponent in the field due to an increase in mass of the mound 
resulting from an increase in accretion rate. Our simulations 
show that for a small change of mound size we have an ap- 
preciable change in maximum magnetic field at the top of 
the mound (Fig. HJ, e.g. between Z c of 55 m and 60 m, the 
maximum field at the top of the mound changes by ~ 36%. 
Since we work in the static limit and do not consider ac- 
cretion rate as a parameter in our simulations, we cannot 
directly probe the luminosity dependence of the spectrum. 
However we can conclude that small changes in height of 
the mound can result in significant changes in the magnetic 
field inferred from the CRSF. 

(vii) Effect of anisotropic continuum : We have also 
carried out runs with mildly anisotropic continuum inten- 
sity profiles : I {6 at) = ^o(l + cos8 a i) and have found no 
appreciable phase dependence of the line energy of a one 
pole CRSF although the percentage modulation of the flux 
was larger. However if continuum is highly anisotropic, e.g. 
I = (sin 20 a i I (29 a i)) 2 then the resultant spectra are found 
to be phase dependent. However such strong beaming may 
not be realistic in the context of HMXB pulsars. 

(viii) Screening effect from overlying column : Ef- 
fects of partial screening of the mound by an atmosphere 
of accreted matter have not been considered in our simu- 
lations. An optically thick blanket of settling plasma can 
screen a fraction of the mound depending on the viewing 
geometry. The effect of screening can be estimated approx- 
imat ely by using the veloc ity profile of the settling flow 
from iBecker fc Wornl (|2007l l : v(z) = 0A9cyG~/Rp~, for a 
neutron star of mass ~ IAMq and radius ~ 10km, mean 
plasma temperature ~ lOkeV and magnetic field ~ 10 12 G. 
For an accretion rate M ~ 10 16 gs _1 and a mound of height 
Z c ~ 70m (maximum allowed Z c for eq. 1 1 1 p we get the 
Thomson optical depth as r ~ 1.57 x 10 _5 £, where £ is the 
line element along the los through the accretion column. For 
a hotspot of radius R p ~ 10 5 cm, we see that regions near the 
axis will be optically thick along a radial line from the axis 
(in local cylindrical geometry of the mound). The degree of 
screening of different parts of the mound will vary with spin 
phase, resulting in phase dependent spectra. 

Thus we conclude from this work that the distortions in 
local magnetic field due to confinement of accreted matter 
at the polar cap has considerable influence on the phase 
resolved spectra from accreting binary neutron star systems. 
Current observations often have poor count statistics and 
involves averaging in phase due to which many of the details 
in the spectra may be lo st. Future missions like ASTR OSAT 
(Agarwal et. al. 2005 , iKoteswara Rao et all |jooj )) with 
higher sensitivity and better energy resolution can help us 
analyse the spectra in greater detail and provide us clues 
to the nature, shape and geometry of the accretion mound. 
A more detailed analysis including dynamic simulations of 
the emitting region and full radiative transfer would help 
us investigate effects like line energy-luminosity correlation, 
anharmonic separation of the CRSF energies, asymmetric 
shape of the CRSF etc. 
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APPENDIX A: SOLVING THE 
GRAD-SHAFRANOV EQUATION 

Al The numerical algorithm 

We express the G-S equation in non-dimensional form by 
scaling the physical parameters with Lo — R p = 1km, po = 
10 6 g cm -3 , Bo = 10 12 G. The flux function is normalised 
to the value at r — R p : u — ip/ip p (ip p = ^BoR p ). We 
perform a variable transformation x — r 2 and solve the G-S 
equation on a (x — y) grid which are the normalised (r 2 — 
z) coordinates. The result is then transformed back to the 
(r—z) grid by spline interpolation. With this transformation 
and scaling, the G-S equation for an adiabatic gas reads as 



if 



<) 2 



Ax d^ + W 2= - Cx{Y °- y) ' 



t dY [u) 
du 



(Al) 



where C = lQitAgLl +cl /Bg, a = ^ and Y (u) is the 
mound height profile expressed in scaled flux coordinate u. 
We adopt the following set of boundary conditions : 

(i) %p = ipd at z = V r. ipd = \Bqt 2 is the initial guess ip 
which is the ip for a uniform field B = Bqz (approximation 
of dipolar field in the polar cap reg ion). This is the line- t ying 
condition of the field at the base ijHameurv et al.lll983T ). 

(ii) ip — ipd at r — V z. 

(iii) ip = ipd set r — » oo V z. Field is undistorted for 
r >> R p . Ideally one should keep the r = r max boundary 
at infinity to fully capture the distortion of the field lines 
resulting from the lateral pressure of the confined plasma. 
To implement this numerically is impractical. In our GS- 
solver the boundary along the radial direction is chosen at 
r m ax = fp where ip = ipd is fixed as the boundary condition. 
The mound height profiles considered for the GS-solver fall 
off with radius vanishing at r — r p . Thus at r — r p bound- 
ary there is very little deviation from initial unloaded field 
configuration. Possible limitations of this are discussed in 
Sec. ELl 

(iv) ip = ipd at 2 — > co. Field is undistorted for z — > oo. 
Setting boundary height (ztop) too close to the maximum 
height of the mound (Z c ) affects the solution and gives in- 
correct result. It is seen that setting z t0 p > 1-5Z C is sufficient 
for the solution to be stable. Further change in boundary 
height (ztop) does not change the solution significantly. 

We adopt a numeri cal sc heme similar to the one fo llowed 
by iMouscho^iasI (|l974T ) and iPavne fc Melatosl (|2004T ). The 
density is determined from the mound height profile and 
eq. 0. A two state iterative scheme using an outer under 
relaxation scheme and an inner Successive Over Relaxation 
(SOR ) method using Chebyshev acceleration (|Press et al.l 
Il993h is adopted to tackle the non-linear source term on the 
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Figure Al. Mean error vs number of iterations of GS-solver till 
fOO iterations. Convergence (eq. I A2 1 was reached after 9 steps. 



RHS. We start with a initial guess solution (ipi n i ~ ip&)- 
SOR loop converges when £ < e S or€in i where £ is the re sidue 
of differenced equation (LHS-RHS) (|Press et al.lll99ot ). £i n i 
is the initial residue at the start of the G-S iterations and 
e S or is the error tolerance for the SOR scheme. The solution 
from the SOR at the i th G-S iteration is then evolved by an 
under-relaxation scheme 



A = 4> l - 



(A2) 



where £ < 1 is the under-relaxation parameter. The conver- 
gence is considered achieved when maximum error at each 
grid point is reduced below threshold 



ip* - -ip*- 



< £GS 



(A3) 



Smaller values of £ gave faster convergence but convergence 
rate did not improve much with £ below an optimum value 
~ 0.01. The convergence rate depends on the model being 
used. Cases for which physical parameters (magnetic field, 
pressure etc) are near regions where the solution does not 
exist (see Sec. I3.3[) . the G-S loop takes longer time to con- 
verge. Otherwise convergence is reached within less than 100 

steps. Fig. lAll shows mean error ( ^ T ~ ) at each step 

N X ,N V 



N X N V 



of iteration for 100 steps. The convergence was reached after 
9 steps in this case. The error tolerance limits were usually 
set at e sor = 10~ 8 and e G s = 10~ 7 for a grid of 1024 x 1024. 



A2 Testing the GS-solver 

Our GS-solver was tested by comparing results with analyt- 
ical solutions of the G-S equation. We solved the Soloviev 
equation 



A 2 tfj = r 2 



(A4) 



whose analytical solution is i/i ana = ~r z ■ The maximum 
absolute difference of the solution from the GS-solver and 
the analytic expression was 



IpB. 



-\ < 9.85 x 10" 



(A5) 



which we consider acceptable. 

We also tested the equili bria by using the MHD code 
PLUTO (|Mignone et al.ll2007l ) and checking if the solution 
is stationary with time. The typical Alfven time scale of 




Figure Bl. The position and orientation of the coordinate sys- 
tem local to the hotspot (X",Y",Z") with respect to the global 
coordinate system (X,Y,Z) defined in the frame of the neutron 
star. The Z axis is the axis of rotation and Z" is along the mag- 
netic pole O' (which is also the centre of the hotspot). The mag- 
netic axis (Z") is at an angle r\ v to the spin axis (Z). The line 
of sight makes angle i with the Z axis and u> (phase) with the Y 
axis. In the frame (X",Y",Z"), (p, <j>, £) are the cylindrical coor- 
dinates where angle <j> is measured with respect to X" axis. X and 
X" axes are in the same direction. 



the accreted mound with a scale length Lo = 1 km, den- 
3 and magnetic field Bo = 10 12 G is 

1-4 



sity po = 10 g cm 

tA = L /Va ~ 3.5 x 10~ 4 s (Va is the Alfven velocity). 



Solution from the GS-solver was put as initial condition in 
PLUTO with a simulation region located inside the density 
mound. At the boundary, values of the physical parameters 
(pressure, magnetic field, density) were kept fixed to the ini- 
tial value obtained from the G-S solution. The simulation 
was run for several Alfven times and the solution was found 
to be stationary. This confirms that the solution obtained 
from our GS-solver represents an equilibrium. 



APPENDIX B: GEOMETRY OF HOTSPOT 
WITH RESPECT TO LINE OF SIGHT (LOS) 

Fig. IB1I shows the relative position of the hotspot with 
respect to the observer and the neutron star. We define 
the frame (X, Y, Z) with origin at the centre of the neu- 
tron star (O) such that X" and X are in the same direction 
and angle between Z,Z" is r\ v . The coordinates of the line 
of sight (los) with respect to the origin O are (i,oj) with ui 
measured from the Y-axis. The unit vector along the los is 
= (sini sincd, sini cos uj, cos i). The hotspot lies on the 
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Figure B2. The relative orientation of the local normal fi;, local 
direction of emission h a of a ray of light and the direction of the 
bent ray n^, which is along the los in the local coordinate system 
of a representative mound. The line along f is the direction of 
the radius vector from the star centre (O) to the mound (see Fig. 
IBlJl , A ray of light along n a makes an angle a with the radius 
vector f which due to gravitational bending is bent towards n^. 
The local emission direction makes an angle 9 a [ with the local 
normal. The vectors r, h a , n^, and h[ all lie in the same plane. 



surface of the star with its centre at O' whose coordinates 
with respect to the O are (R s ,rj p ) (see Fig. IB1[I . Z" being the 
normal through O'. The G-S computation is done in a cylin- 
drical coordinate system (p, </>,£) in the frame (X",Y",Z") 
where angle <p is defined with respect to X". Coordinates of 
a point in (X",Y",Z") frame is (pcos <p, psin </>, £). 

To express the coordinates of a point on the hotpost in 
the (X,Y,Z) frame, we first rotate the (X",Y",Z") about 
X" by angle r/ p to (X', Y',Z') such that Z' is parallel to Z, 
and then shift the origin from O' to O by distance R 3 along 
the radial line joining OO'. Thus we get the coordinates 
of a point r on the hotspot as r = {pcos (j>, pcos r\ v sin</> + 
(C + Rs ) sin r/p , (£ + R s ) cos r/ p — p sin rj p sin <p} , where r = 



(p 2 + (£ + Rs) 2 ) 1,/2 is the radial distance of (x,y,z) from 
O. 

The angle between the los (n^,) and the radius vector to 
a point on the hotspot is cosi/> = ny, • r/r. The ray coming 
towards the observer along a los is deviated from its local di- 
rection of emission on the neutron star surface (n a , see Fig. 
IB2|) due to gravitational bending of light. The local emis- 
sion angle (a) and the angle between the los and the radius 
vector f can b e related by the approximate formula given by 
Belob orodov (|Beloborodovl |2002| ; iPoutanen fc Beloborodovl 
120061 ) : 



+ (1 — u) cos 



(Bl) 



where u — r 3 /r, r s being the Schwarzschild radius. Since 
the unit vectors n a , n,/, and f lie in the same plane, we can 
write n Q x (n a x f ) = Cn a x (ny, x f), C being a constant, 
as the two vectors lie in the same direction differing only in 
magnitude. Using h a ■ r = cos a, the constant C is evaluated 
to be C = sing/ sirup, for a ^ 0. So n a ca n be related to 
ny, (as in lPoutanen Sz Gierliriskil (|2003h ) 



sin(i/j — a) a 
- r 



sin a „ 



(B2) 

sin ip sin tp 

For a = (which corresponds to ip = from eq. IB1|) . 
we have r\ a = = f . Numerical computation of sin(^/j — 
a)/ sim/) and sin a/ sintp leads to large errors in the limit 
a) —¥ 0. We expand eq. ( IB1[) for small a and ip : a = 



— u)ip. Using this, we get the following approximate 
relations 

sin(?/> — a) 



sinip 
sin a 



sin?/) 



u + cos ipVT^u(VT~ 7 u- 1) (B3) 
y/T^L (B4) 



The errors in the approximate form of eq. (|B2I) for small 
(ip,a) are less than 3% for eq. (|B3|) and 1.5% for eq. (|B4|) 
for V < 25°. 

The angle between the local direction of emission (n a ) 
and the local magnetic field (b) is required to determine 
the angular dependence of the width and relative intensity 
of the cyclotron resonance scattering features (CRSF) (Sec. 
14. ip . From eq. (|B2|) and unit vector along local magnetic 
field b = b p p + we get the angle a t between n Q and b 
as : 

cos 6 a b = n a • b = Ci cos 6 r b + Ci cos Q^b (B5) 

where cos Q r b = f • b = (pb p + (^ + R s )b^)/r, and C\ and C2 
are the coefficients of f and respectively from eq. (|B2|) . 
eq. (|B3|I and eq. (|B4|I whichever is applicable. 



cos 6^b = b 

= bp [sin i sin ui cos <p + (sin i cos to cos r\ v — cos i sin r\ v ) sin p] 
+ [sin i cos sin r\ v + cos j cos r]p] 

The angle between n a and the local normal (n;) is 
required to evaluate the flux from a local area element 
(eq. I18[) for a given intensity profile. To find the unit 
vector along the local normal to the mound (rii) we first 
fit the top profile of the mound (£top = f(p) □) with a 
polynomial function of p (the order of the polynomial is 
chosen to keep errors to less than 5%). From the slope 
(m s — d^top/dp) we find the normal to the mound in the 
local frame (X",Y",Z") and transform it to the (X, Y, Z) 
coordinate to get n; = {— sin 9 S cos <p, — sin 9 S cos r/ p sin <j> + 
cos 6 a sin T] p , cos a cos rj p + sin 9 S sin r/p sin <f>} , where sin 9 S = 



and cos 9 S = 



between n a and n; as : 



Thus we get the angle 9 a 



cos a i — n a ■ hi — Ci cos 9 r i + C2 cos 9 



(B6) 



where again Ci and C2 are the coefficients of f and ny, re- 
spectively from eq. (|B2p . eq. (|B3p and eq. ()B4|I whichever 
is applicable. cos# r ; = r ■ n; = (£ + R s — m s p) / (ry/1 + m'j) 
and 

cos Qt/,1 = £i,p ■ ni 
m 



[(cos i sin r\p — sin i cos to cos r; p ) sin < 



vTTm 



sin i sin oj cos + 



[sin i cos w sin 77^ + cos i cos rj p \ 
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